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We present gravitational waveforms for the last orbits and merger of black-hole-binary (BBH) systems along 
two branches of the BBH parameter space: equal-mass binaries with equal non-precessing spins, and nonspin- 
ning unequal-mass binaries. The waveforms are calculated from numerical solutions of Einstein's equations 
for black-hole binaries that complete between six and ten orbits before merger Along the equal-mass spinning 
branch, the spin parameter of each BH is Xi = Si/MJ £ [—0.85,0.85], and along the unequal-mass branch the 
mass ratio is (/ = M2/M1 e [1,4]. We discuss the construction of low-eccentricity puncture initial data for these 
cases, the properties of the final merged BH, and compare the last 8-10 GW cycles up to Mco = 0.1 with the 
phase and amplitude predicted by standard post-Newtonian (PN) approximants. As in previous studies, we find 
that the phase from the 3.5PN TaylorT4 approximant is most accurate for nonspinning binaries. For equal-mass 
spinning binaries the 3.5PN TaylorTl approximant (including spin terms up to only 2.5PN order) gives the most 
robust performance, but it is possible to treat TaylorT4 in such a way that it gives the best accuracy for spins 
Xi > —0.75. When high-order amplitude corrections are included, the PN amplitude of the {(. = 2,m = ±2) 
modes is larger than the NR amplitude by between 2-4%. 
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I. INTRODUCTION 

One of the most urgent goals of numerical relativity is to 
produce simulations that will aid the detection of gravitational 
waves (GWs) from black-hole -binary mergers. The current 
first generation of ground-based interferometric GW detectors 
is about to be upgraded, and the second-generation Advanced 
LIGO and Virgo detectors are expected to come online around 
2014 |1-4|. Once operational, current event-rate calculations 
predict that they may observe multiple GW signals in one 
month of design-sensitivity operation [l]. Some of these sig- 
nals will be from the inspiral and merger of two black holes, 
and to find them in the detector data GW astronomers will use 
matched filtering techniques, for which they require large col- 
lections of accurate theoretical waveforms (templates) of the 
physical signal. 

The GW signal from the last orbits and merger of black- 
hole-binary systems can only be calculated in full general rela- 
tivity using numerical solutions of Einstein's equations. Since 
such simulations became possible in 2005 |5-7|, they have 
been used to explore larger regions of the black-hole-binary 
parameter space (which is parametrized by the mass ratio of 
the binary, the spin vector of each black hole, and the binary's 
eccentricity), with increasing levels of accuracy and covering 
increasing numbers of GW cycles before merger (SI |9J . In 
addition to use in producing analytic waveform models for 
the construction of GW search template banks, which we will 
discuss fuither in this paper, NR waveforms have also been 
useful in GW detection efforts as part of the NINJA project to 
test a battery of current GW search pipelines ifToHTTI . 

In this paper we present simulations that cover between 
six and ten binary orbits before merger of configurations in 
two important sub-families of the binary parameter space: 
unequal-mass binaries in which the black holes are not spin- 



ning, and equal-mass binaries where the black holes have 
equal spins either aligned or anti-aligned with the binary's or- 
bital angular momentum. 

Following a brief summary of numerical methods in Sec. In] 
in Sec. Ill we extend the method we developed in [12 1 to pro- 
duce low-eccentricity parameters for spinning binaries. This 
method is based on integrating the PN equations of motion 
from a separation where quasi-circular parameters are suffi- 
ciently accurate, up to the binary separation where we wish to 
begin a full numerical simulation, and then using the momenta 
at that separation from the PN integration as the initial mo- 
menta of the full numerical simulation. We now incorporate 
the highest-order known spin contributions, but find that these 
are still not accurate enough, and develop a method to further 
refine the PN predictions. This allows us to produce simula- 
tions of unequal-mass nonspinning and equal-mass spinning 
binaries with eccentricities of e < 0.004. 

In the Samurai study |fT3l it was shown that cuiTent numer- 
ical simulations for the equal-mass nonspinning case are well 
within the accuracy requirements for detection with ground- 
based experiments. That study also showed that the agreement 
of numerical results between different codes was consistent 
with the error estimates of each code — and so a complete er- 
ror analysis of a set of numerical simulations can confidently 
be considered as providing the uncertainty in those waveforms 
with respect to the true physical waveforms. In Sec. [V] we 
study the eiTors in our unequal-mass and equal-mass-spinning 
waveforms, and conclude that these waveforms are also well 
within the accuracy requirements for GW detection. 

We also estimate the phase accuracy of our simulations, us- 
ing a number of different methods. The phase eiTor accumu- 
lates quickly during the inspiral and even faster during the 
merger, and small eiTors at a given separation or frequency 
are amplified during the further evolution. On the other hand. 



the absolute value of the GW phase is not directly observable, 
and phases from different simulations can be aligned in differ- 
ent ways, e.g. between two suitably chosen fixed frequencies 
during the evolution. Correspondingly, the estimated phase 
errors show a dramatic dependence on such alignment effects, 
and it is therefore useful to phrase error estimates in a number 
of different ways. An example relevant to a comparison with 
PN results would be the time-domain phase error over a given 
set of GW cycles, while for GW detection we might be more 
interested in the mismatch-error in the waveform at a given 
binary mass with respect to a given detector 

Having described the production of our numerical wave- 
forms, and established their accuracy, we summarize in 
Sec. |VT]the physical properties of the configurations we have 
studied: the mass and spin of the final merged black hole, the 
recoil (in the unequal-mass cases), and the radiated energy 
and its distribution among the dominant and sub-dominant 
harmonics. 

One of the motivations for producing black-hole-binary 
simulations is to use the resulting waveforms for the con- 
struction of analytic waveform models. All such models are 
based in some way on the expectation that PN approxima- 
tions will be sufficiently accurate up until a few orbits before 
merger, and so PN (or effective-one-body) results can be used 
to model the early inspiral, and numerical results can be used 
to calibrate a model for the merger. Such a procedure first 
requires a measurement of the accuracy of PN results as the 



binary approaches merger. We do this in Sec. VII where we 



compare the phase and amplitude of our NR results with the 
corresponding PN predictions over the 8-10 cycles prior to the 
point where the GW frequency reaches Mco = 0.1, about 1.5 
orbits before merger. This extends our previous studies of the 
equal-mass nonspinning binary |14, 15 1 and equal-mass bina- 
ries in the orbital-hangup configuration llT6l . 



II. NUMERICAL METHODS 

We performed numerical simulations with the BAM code 
ifTTl [TSl . The code starts with black-hole-binary puncture 
initial data flW, "201 generated using a pseudo-spectral ellip- 
tic solver 1 21 1, and evolves them with the ;i^-variant of the 
moving -puncture El?] version of the BSSN Il22ll23l formu- 
lation of the 3h-1 Einstein evolution equations. Spatial finite- 
difference derivatives are sixth-order accurate in the bulk ifTSl . 
Kreiss-Oliger dissipation terms converge at fifth order, and a 
fourth-order Runge-Kutta algorithm is used for time evolu- 
tion. The gravitational waves emitted by the binary are calcu- 
lated from the Newman-Penrose scalar ^4, and the details of 
our implementation of this procedure are given in fTTl. 

In each simulation, the black-hole punctures are initially a 
coordinate distance D apart, and are placed on the y-axis at 
yi — —qD/{\ +q) and >'2 ~ D/{\ +q), where q — M2/M1 
is the ratio of the black hole masses in the binary, and we 
always choose Mi < M2. The masses M; are estimated from 
the Arnowitt-Deser-Misner (ADM) mass at each puncture, ac- 
cording to the method described in |19|; we discuss this esti- 
mate of the black-hole masses and its subtleties for spinning 



black holes in more detail in Appendix IA] The Bowen-York 
punctures are given momenta p^ = ^pt tangential to their sep- 
aration vector, and py = zLp^ towards each other. The latter 
momentum component accounts for the (initially small) ra- 
dial motion of the black holes as they spiral together. One 
essential question in setting up our simulations is the deter- 
mination of the parameters {p,,p,-) that lead to non-eccentric 
quasi-circular inspiral. We will discuss our procedure to gen- 



erate low-eccentricity parameters in Sec. Ill 



The grid setup is similar to what we have used in fTTl, and 
using the notation introduced there, the simulations discussed 
in this paper all use a configuration of the form Xmti=2 [h 'xN : 
h X 2N : 6]. This indicates that the simulation used the 
X variant of the moving-puncture method, /i nested mesh- 
refinement boxes with a base value of A^^ points surround each 
black hole, and I2 nested boxes with (2A^)^ points surround the 
entire system, and there are six mesh-refinement buffer points. 
The Tj parameter in the BSSN system is Mr\ — 2. The choices 
of A^, l\, h and the resolutions are given in Tab. |l] The res- 
olution around the puncture is denoted by Mi //!„„„, which is 
the resolution with respect to the smallest black hole, M\ . The 
puncture of the second black hole will have the same numeri- 
cal resolution, but if the black hole is bigger, M2 > Mi, then it 
will effectively be better resolved. 

The one exception to this setup is a second convergence se- 
ries for mass ratio q — A- (see the last row in Tab.ll|i. These sim- 
ulations use a grid configuration in which the effective finest 
resolution is the same for both black holes. This is achieved 
by putting different numbers of refinement boxes around each 
puncture. As Mi is four times smaller than M2, we use two 
more boxes (the resolution doubles from box to box) around 
the smaller black hole than we do for the larger one. 

Far from the sources, the meaningful length scale is the total 
mass of the binary, M = Mi +M2, and so the resolution on the 
coarsest level is given by h^ax/M. 

In our previous study of Xi > cases [16| we found that 
extra resolution was required around the punctures when the 
black holes have high spin, \Xi\ > 0.75. In the newer Xi < 
simulations we use high resolution around the puncture in 
all cases. Note also that in only two cases (q = 2,3), is the 
outer boundary causally disconnected from the physical sys- 
tem for the entire length of the simulation. This can be seen 
by comparing the time when the GW signal reaches its peak 



amplitude in Tab. IV and the location of the outer boundary 
in Tab.!] 



III. SPECIFICATION OF LOW-ECCENTRICITY INITIAL 
PARAMETERS 

Puncture initial data typically consist of the analytic 
Bowen-York solution to the momentum constraint f20l (which 
allows for the construction of multiple boosted, spinning black 
holes), and a numerical solution of the Hamiltonian constraint 
in puncture form 1 19 1. To produce the data for the simulations 
that we discuss in this paper we used the single-domain spec- 
tral elliptic solver described in |21 1. In this approach the black 
holes' angular and linear momenta may be directly specified. 



TABLE I. Summary of grid setup for numerical simulations. The grid parameters follow the notation introduced in 1171 : see text. M\/h,„i„ 
denotes the resolution on the finest level with respect to the smallest black hole, while h,„ax/M is the resolution on the coarsest level with 
respect to the total mass, M = Mi +M2. The outer boundary of the computational domain is at Xij„ax/M, where x, = {x,y,z}- In general /; 
indicates the number of moving refinement levels around each puncture, and I2 the number of large refinement levels that encompass both 
punctures. The one exception is the second q = 4 series, which uses three refinement levels around the puncture of the large black hole, and 
five around the other. 



Configuration 


A' 


{luh) 


Ml /hmin 


hmax/M 


Xi,nmx/M 


Equal-mass simulations 


Xi = -0.85 


72,80,88 


(6,5) 


48,53.3,58.67 


10.67,9.6,8.73 


114 


Xi = -0.75 


80,88,96 


(6,5) 


53.3,58.67,64 


9.6,8.73,8.0 


774 


Xi = -0.50 


64,72,80 


(6,5) 


42.67,48,53.3 


12.0,10.67,9.6 


774 


Xi = -0.25 


64,72,80 


(6,5) 


42.67,48,53.3 


12.0,10.67,9.6 


774 


Xi^O 


64,72,80 


(5,5) 


21.3,24,26.7 


12.0,10.67,9.6 


774 


Xi = 0.25 


80,88,96 


(5,5) 


26.7,29.33,32 


9.6,8.73,8.0 


774 


Xi = 0.50 


80,88,96 


(5,5) 


26.7,29.33,32 


9.6,8.73,8.0 


774 


Xi = 0.75 


64,72,80 


(6,5) 


42.67,48,53.3 


12.0,10.67,9.6 


774 


Xi = 0.85 


64,72,80 


(6,5) 


42.67,48,53.3 


12.0,10.67,9.6 


774 


Unequal-mass 


simulations 










q = 2 


70,80,88 


(5,7) 


23.3,26.67,29.33 


29.26,25.6,23.27 


2063 


q = 3 


70,80,88 


(5,7) 


23.3,26.67,29.33 


21.94,19.2,17.45 


1547 


q = 4(a) 


70,80,88 


(5,7) 


23.3,26.67,29.33 


17.55,15.36,13.96 


1237 


g = 4(b) 


80,88,96 


(3/5,7) 


26.67,29.33,32.0 


15.36,13.96,12.8 


1237 



while their masses are specified indirectly. Some subtle issues 
related to the estimation of black-hole masses are discussed in 
Appendix [A| 

We wish to simulate black holes following non-eccentric in- 
spiral. Eccentricity cannot be so easily defined in full general 
relativity as in Newtonian theory, since radiation reaction pre- 
cludes the existence of circular orbits, and gauge effects mean 
that any definition based on the coordinate motion of the black 
holes (or the punctures) will not be unique. Having said that, 
all definitions of eccentricity based on the gauge-invariant 
gravitational-wave signal (see |24| for a thorough discussion 
of the choices available) should agree on zero-eccentricity in- 
spiral, and we have found in previous work |12j that defini- 
tions based on the coordinate motion are acceptable at the 
level of accuracy that we are interested in, which is eccen- 
tricities on the order of e ^ 10^^. For this work we estimate 
the eccentricity from the orbital frequency of the puncture mo- 
tion by the maximum of 6^, = (a)(f) — a)c(f))/(2a)f(f))' where 
(i)c{t) is an estimate of the non-eccentric frequency based on a 
curve fit through the numerical data lfT2l . 

To produce low-eccentricity inspiral we need to know the 
appropriate initial momenta to give the black holes. In punc- 
ture simulations the most effective way to do this seems to 
be to estimate the momenta using PN theory. The simplest 
approach is to consider only a conservative PN model (i.e., 
without radiation reaction, and therefore without inspiral), and 
to calculate the momenta consistent with circular orbits at a 
given coordinate separation. We will refer to these as quasi- 
circular (QC) parameters. We have used QC parameters in the 
past ifTTl I25I |26 I for binaries no more than three orbits away 



from merger, at which point small eccentricities are hard to 
detect. For binaries that undergo five or more orbits before 
merger, it becomes clear that the QC momenta result in no- 
ticeable eccentricities. 

One way to improve the parameters is to use a PN approach 
that includes the effects of radiation reaction. A straightfor- 
ward method to do this is by time integration of the PN equa- 
tions of motion. One begins with QC parameters for a bi- 
nary with a large separation (D > 30M), and integrates the 
PN equations of motion for two point particles until they have 
reached the separation at which we wish to start a full numer- 
ical simulation. At that point we read off the parameters from 
the PN calculation and use those in our black-hole evolution 
code. We will refer to such parameters as PN inspiral (PN) 
parameters. In [121 we demonstrated that, in the equal-mass 
nonspinning case, PN inspiral parameters using a 3PN accu- 
rate Hamiltonian 127^291 (see also Il30ti32l ) and 3.5PN accu- 
rate radiation flux ||33ti35]| lead to inspirals with eccentricities 
e ~ 0.002. Similar results were also obtained for unequal- 
mass nonspinning binaries up to mass ratio q = M2/M1 = 4; 
these were first used in fSSl, and are discussed in more detail 
here. PN inspiral parameters were also used successfully to 
produce low-eccentricity simulations of one precessing-spin 
configuration in [37|. 

When we extended our studies to spinning binaries in lfT6l 
we found that the same procedure did not work so well. At 
that time the PN equations included only leading-order (LO) 
contributions to the spin-orbit and spin-spin Hamiltonians 
[38-43 1, and spin-induced radiation flux terms as described 
in 1,44] . see also 1.41., 43J . Note that there can be ambiguities 



TABLE II. Choices of the initial momenta and resulting eccentricity 
for the equal-mass Xi = —0.5 case. We find that the quasi-circular 
(QC) parameters yield an eccentricity an order of magnitude larger 
than we desire, while the leading-order (LO) PN inspiral parameters 
are twice as bad. Incorporating NLO spin effects dramatically re- 
duces the eccentricity to e ~ 0.004, and a further iteration based on 
PN predictions reduces it by a further 25%. 



Parameters 


Px 


Py(xl0-4) 


e 


QC 


0.08469 





0.015 


PN (LO) 


0.08612 


-5.824 


0.03 


PN (NLO) 


0.08500 


-5.250 


0.004 


PN+ 


0.08512 


-5.258 


0.003 


PN_ 


0.08487 


-5.242 


0.008 



in the literature assigning PN orders to the spin-terms in the 
Hamiltonian. Following |45| we assign 1.5PN order to the 
LO spin-orbit term in the Hamiltonian and 2PN order to the 
leading-order spin-spin term. 

We find that the corresponding LO PN inspiral parameters 
lead to inspirals with eccentricities e ^ 0.009 for spins parallel 
to the orbital angular momentum (and Xi = Sj/Mf — 0.5), and 
e ~ 0.03 for the corresponding case with spins anti-parallel 
to the orbital angular momentum. In fact, in many cases 
much lower eccentricities were achieved with the supposedly 
cruder QC parameters, and these were used for the final re- 
sults in lfT6ll : in spinning cases the QC parameters include only 
leading-order spin effects ll26ll4TI . and we have not examined 
the effects of using NLO QC parameters. For the anti -parallel- 
spin cases, both the QC and LO PN inspiral parameters lead 
to eccentricities that were too high to be seriously considered 
as "quasi-circular inspiral". 

For this work we have incorporated recent results ll45H47l . 
and include next-to-leading order (NLO) spin terms in the PN 
equations of motion. We have also included the flux contri- 
bution due to the energy flowing in to the black holes, which 
appears at the relative 2.5PN order, as derived in Ref. |48|. 

The improvement when using NLO PN inspiral parameters 
is dramatic, as shown in Table. [Ill ^^ ^^^ ^^^^ ^^ '^^e ;i^, = — 0.5 
case, the QC parameters lead to an eccentricity of e w 0.015. 
The LO PN inspiral parameters give even higher eccentricity, 
e « 0.03. When NLO spin terms are included, however, the 
eccentricity drops to e « 0.004, i.e., a reduction by almost an 
order of magnitude. 

However, this eccentricity is still twice what we can achieve 
in the nonspinning case. One approach to reduce further the 
eccentricity would be to employ an iterative procedure like 
that used for excision data in |49|. We do not attempt this 
procedure, for the following reasons. The excision data con- 
sidered in 1 49 1 are adapted to the gauge that will be used for 
their subsequent evolution, and in particular already possess 
the coordinate motion consistent with their motion along a 
quasi-circular inspiral. Puncture initial data, by contrast, start 
out with no coordinate motion. After the simulation begins, 
the puncture wormholes evolve into puncture trumpets ll50l - 
|32 1, and acquire some coordinate motion which, after roughly 
one orbit, corresponds to the motion consistent with quasi- 



circular inspiral fTZJ . Measuring the orbital eccentricity in or- 
der to apply an iteration procedure therefore requires perform- 
ing puncture simulations well beyond one orbit, and even then 
the reduction in eccentricity typically converges very slowly. 
In addition, it is only really practical to estimate the eccentric- 
ity using the puncture coordinate motion, which, as should be 
clear from the preceding discussion, is purely a gauge effect 
and so may not be able to be used to reduce the physical ec- 
centricity to an arbitrarily low level. It turns out that a similar 
effect occurs with excision data, as recently reported in fl^. 
Note, however, that the results in ll24l arise from a different 
gauge condition to that used in moving-puncture simulations, 
and any conclusions they draw about the correspondence (or 
lack thereof) between the coordinate and physical motions of 
the black holes may not apply to moving-puncture results. We 
will consider this further in future work |53 1. 

For these reasons, we have attempted an alternative ap- 
proach, based on extracting further information from PN the- 
ory. We will illustrate our approach using the same case 
that we have discussed above. Here the black-hole spins are 
\Xi\ — 0.5 and directed anti-parallel to the binary's orbital an- 
gular momentum. The initial coordinate separation of the 
punctures is D = 12. 5M. A first simulation is performed us- 
ing NLO PN inspiral parameters. For reference, these are, as 
given in Table [II] pt/M — 0.0850 for the tangential momenta, 
and pr/M = —5.250 x lO^'* for the radial momenta. 

We then calculate the eccentricity of the first orbits in this 
simulation. As we stated above, we find e sa 0.004. To re- 
duce the eccentricity further, we now return to PN theory. If 
we solve the PN equations of motion starting at D = 12.5M, 
using {pr,pt)lM = (0.085,-5.25 x lO""*), we wiU of course 
recover the same non-eccentric PN inspiral as when we first 
calculated these parameters. We now ask the question, "How 
much would these parameters have to vary, in order to pro- 
duce the eccentricity we saw in our NR simulation?" We as- 
sume that a variation in {pr,Pr) that produces an eccentricity 
of e = 0.004 in the PN integration, starting at the same sep- 
aration as the NR simulation, will give us approximately the 
correct magnitude of momenta variation to remove the eccen- 
tricity in the NR simulation. Common simplifying character- 
istics of the situations we have in mind are that the variation 
in separation due to eccentricity, during half an orbit, say, is 
smaller than the variation due to the inspiral; and that the tan- 
gential momentum is much larger than the radial momentum. 

What we do now is simply adjust both the radial and tan- 
gential momenta by some factor k until the PN evolution pro- 
duces an eccentricity of e w 0.004. For this case we find that 
k K, 1.0015, i.e., by 0.15%. We then assume that this is close 
to the error in the parameters that we have used in our NR 
simulation, and modify those also by 0.15%. We perform 
two additional simulations, one in which the momenta are in- 
creased by 0. 15%, and one in which the momenta are reduced 
by 0.15%. The two enhanced PN-parameter choices are de- 
noted by "PN±" in Table[II| 

We find, in this case, that the eccentricity is reduced when 
the momenta were reduced, and we achieve e « 0.003. This 
is very close to the eccentricity we achieved in the equal-mass 
nonspinning case, and we consider that acceptable. We also 



experimented with repeating this iteration procedure, but it 
did not noticeably reduce the eccentricity; at this level it may 
be possible to more cleverly modify separately the radial and 
tangential momenta, but it may also turn out that further re- 
finement in the puncture-motion eccentricity will not improve 
the true physical eccentricity, as indicated by the GW signal. 
At some level of accuracy it will also be necessary to adjust 
the radial and tangential momenta by different factors. An- 
other shortcoming of the method we have used here is that we 
do not explicitly use phasing information when adjusting the 
eccentricity of the orbit, e.g. to determine whether momenta 
should be increased or decreased without having to perform 
two further simulations. On reason why this is difficult, is be- 
cause in the initial gauge of the simulation the punctures are 
stationary on the numerical grid, and it takes ^ one orbit for 
their motion to asymptote to a trajectory consistent with the 
physical motion of the black holes. We have found the results 
to be acceptable for all cases we have considered here. Further 
work on improving this procedure is underway, and prelimi- 
nary results on an improved method that also uses phasing 
information have been presented recently f5A\ . 

The procedure we have described here was performed on all 
of the anti-hangup cases, and the final parameters are given in 
Table |IV] For these cases we also indicate the eccentricities 
that were achieved from the raw PN-inspiral parameters. Also 
shown are the parameters for the nonspinning case presented 
in I.14J , the hang-up cases described in I.16J , and a set of non- 
spinning unequal-mass simulations. In the unequal-mass sim- 
ulations, the eccentricity was found to be sufficiently low with 
the raw PN inspiral parameters, and no further modifications 
were made. This also suggests that while we expect the PN 
approximation to deteriorate for larger mass ratios, this dete- 
rioration is not large at ^ = 4. 

However, a second series of ^ = 4 simulations was also per- 
formed, using FOB parameters as described in |55|. Whereas 
PN inspiral parameters lead to an eccentricity of e « 0.0038, 
the FOB parameters lead to a lower eccentricity of only e w 
0.003. This appears to be consistent with the expectation that 
FOB methods retain their accuracy at higher mass ratios bet- 
ter than PN methods, although we note that the uncertainty 
in the eccentricity calculation is 5 x 10^"*, and so the two 
values agree within uncertainty. Results from much higher- 
mass-ratio simulations are necessary to definitively compare 
the performance of FOB and PN parameters. 



IV. SUMMARY OF NUMERICAL SIMULATIONS 

In this section we summarize the two sets of configurations 
that we studied. 

The first comprised equal-mass (^ = 1) binaries, with equal 
spins directed either parallel {Xi > 0) or anti-parallel (Xi < 0) 
to the orbital angular momentum of the binary. The spins con- 
sidered were \Xi\ = {0,0.25,0.5,0.75,0.85}. When the spins 
are (anti-)parallel to the orbital angular momentum, they will 
not precess, making such cases a relatively simple sub-family 
of the total black-hole-binary parameter space. And when 
q = I, the system possesses enough symmetry that the sim- 



ulation can be performed on only a quarter of the full physical 
domain, z > 0,y > 0; this symmetry is also reflected in the fact 
that in these configurations the center-of-mass of the system 
does not move, and the final black hole does not experience 
any recoil. 

In addition, the choice of equal spins also yields an impor- 
tant sub-family of configurations: we found in [561 that it is 
possible to rather accurately model any non-precessing bina- 
ries with unequal spins using essentially only a mass-weighted 
sum of the spins of the binary, motivated by numerical ev- 
idence from |26,!52l|58l, and PN theory |59|. Therefore for 
the purposes of producing waveform models for GW detection 
with current ground-based detectors, it is sufficient to simulate 
only binaries where the black holes have equal spins. Recall 
also that the use of Bowen- York-puncture data limits us to 
black holes with spins \Xi\ < 0.92 |60U63l, and that since er- 
rors due to the presence of junk radiation increase with higher 
spins, the maximum spin that we treat in these simulations is 
lzd-0.85. 

The second set of configurations we consider is nonspin- 
ning binaries with unequal masses, q = {2,3,4}. Now the 
symmetry of the system is reduced, and the simulations re- 
quire half of the physical domain, z > 0. The center-of-mass 
of the system can move, and does due to the nature of the 
asymmetry of the GW emission from these systems, and the 
final black hole "recoils" (or is "kicked") relative to the orig- 
inal center-of-mass of the binary. We will discuss further the 



recoil in our unequal-mass simulations in Sec. VI D 



Of the simulations discussed in this paper, a small sub- 
set were first presented elsewhere. The equal-mass nonspin- 
ning simulations were described in detail in lfT4l . and the 
\Xi\ > cases in |16|. The remaining simulations have not 
yet been published, although they have all be used as part 
of other studies. The q — 2 simulations were used to study 
parameter-estimation accuracy for LISA ["641. In addition, all 
of these simulations have been used to build phenomenologi- 
cal waveform families. Some of the nonspinning-binary data 
were used in ['36' "65 1 to produce nonspinning phenomeno- 
logical waveforms, and in |66| to calibrate an effective-one- 
body (FOB) model. It should be noted, however, that in those 
works less-accurate simulations of the higher-mass-ratio cases 
were used, and in particular the accuracy of the q = A- data 
used in |66| were not sufficiently accurate to conclusively 
test the physical fidelity of the FOB model. Also, all of the 
waveforms presented here were used to produce the first non- 
precessing-spin phenomenological model presented in |56|, 
and the follow-up study in [67.1 . In fact, extra simulations for 
unequal-mass spinning binaries were also necessary for that 
work, but we will not consider those here. 

The methods we used to estimate the initial momenta for 
quasi-circular inspiral were described in Sec. Ill and fall 
into three classes: quasi-circular (QC), PN inspiral (includ- 
ing next-to-leading order [NLO] spin terms when necessary), 
and enhanced PN± inspiral, as described in Sec. 



Ill For an 



extra ^ = 4 series, we also used parameters based on an FOB 
model. QC parameters were used to produce the older orbital 
hangup Xi > simulations, which were originally presented 
in [.16J ; for these cases the eccentricity is the highest, around 



e ~ 0.006. PN inspiral parameters were used for the unequal- 
mass simulations, where they were found to yield acceptably 
low eccentricities of e < 0.004. Finally, PN± -inspiral param- 
eters were used for the anti-hangup cases j, < 0, and for all 
cases the eccentricity is e < 0.003. The BOB parameters that 
were used for the second q~'\ series lead to an eccentricity 
ofesa 0.003. 



V. ACCURACY OF NUMERICAL SIMULATIONS 



servative error estimates are within the error bounds required 
for this work and many GW-astronomy applications. 

We will illustrate these points with two representative 
cases: one that shows clean sixth-order convergence, and one 
that does not. We focus first on the GW phase, because it is 
the phase error that dominates the mismatch calculations that 
are at the heart of the matched-filtering technique used in GW 
searches, and because the error in the amplitude is dominated 
not by numerical resolution but by the radius 7?ex of the GW 
extraction; we will discuss this further in Sec. |VC| 



In this section we will estimate the errors in our numerical 
results. These will give us some indication of both the physi- 
cal accuracy of our waveforms, and the applications for which 
they can confidently be used. 

The waveform can be decomposed into phase and ampli- 
tude functions, (j)(',„{t) and Af„,(f) respectively, for each spher- 
ical harmonic mode {i,m). The individual harmonics of the 
Newman-Penrose scalar ^4 can now be written as 



^ex'P«m(0=^ft«(f)e" 



(0 



(1) 



where i?ex is the coordinate radius of the wave-extraction 
sphere. The frequency of the signal is given by C0e„, = 
d(j)(„,/dt. The results in this paper use GW signals extracted 
at 7?ex = 90M, unless otherwise stated. 

We will focus on the {£ = 2,m = 2) mode. The numerical 
error in the functions A22{t) and 022(0 is estimated by means 
of a convergence test: we perform simulations at three (or 
more) numerical resolutions, and verify that differences be- 
tween successive resolutions decrease at a rate consistent with 
the expected convergence properties of the numerical code. 
Throughout the remainder of this paper, the functions A{t), 
(j){t) and (o{t) will refer to the {£ = 2,m = 2) quantities un- 
less otherwise stated. If a clear convergence rate is observed, 
then it is also possible to use Richardson extrapolation to re- 
move the error term at the next order, and produce a yet more 
accurate estimate of the true result, and to also calculate an 
uncertainty estimate. This procedure was carried out for our 
equal-mass nonspinning data in lfT4l . 

As described in Sec. Ill] we have based our grid configu- 
rations on that used in our work on equal-mass nonspinning 
binaries llT4l . For that configuration, the defining number of 
points in the three convergence series simulations (see Tab. |l]i 
was A^ = {64,72,80}. This was sufficient to achieve reason- 
ably clean sixth-order convergence, which had been identified 
in |,18J as the dominant order of finite-difference error in our 
code. For configurations with higher mass ratios and non-zero 
spins, greater numerical resolution is required. In many of the 
new simulations the number of grid-points has therefore been 
increased. In some cases the extra resolution was sufficient 
to again achieve sixth-order convergence, but in others it was 
not. For these latter cases, although we can be confident that 
the simulations are converging towards the continuum solu- 
tion, we are unable to use Richardson extrapolation to esti- 
mate the uncertainties, and must provide much more conser- 
vative error estimates. Although we could perform more sim- 
ulations at yet higher resolution, we find that even these con- 



A. GW phase 

In studying the GW phase, we have two aims: (1) to show 
that our results are converging to the continuum solution as a 
function of numerical resolution, and to provide error bounds 
on the GW phase, and (2) to illustrate the ambiguities inherent 
in estimating the phase error The ambiguity in estimating 
phase errors is already well known, but we discuss it further 
here to make the point that although the accuracy of the GW 
phase is important, any one method of estimating the "phase 
error" may tell us little about the waveform's accuracy for a 
given application. 

We consider three measures of the phase error The first is 
the total accumulated phase error over the length of the simu- 
lation. This is the most natural quantity to study in a conver- 
gence test. The second is the accumulated phase error over the 
ten cycles up to a GW frequency of Mco = 0.1. We will need 
this error estimate to justify the comparison with PN results in 



Sec. VII but we will also see that this quantity is problematic 



when used for a convergence test, as is any realignment of the 
GW phase. None of these estimates have a natural interpre- 
tation for applications in GW searches. There the mismatch 
error is the appropriate measure of the waveform's accuracy, 
and we will consider this in Sec. IV CI 

For our convergence analysis we consider in detail two 
cases, {q — l,Xi = 0-5) and the second {q = 4,Xi = 0) con- 
vergence series. These are indicative of the general features 
of all of the cases we have studied. 

Standard convergence plots of the GW phase are shown in 
Fig. [T] The initial phases agree at f = 0, and the plot shows 
the subsequent evolution of the phase disagreement between 
simulations at different resolutions. In the Xi = 0.5 case (left 
panel), which uses relatively high resolution for a moderate 
spin value, we see reasonably clean 6th-order convergence. 
The q = 4 case, however, even though the numerical reso- 
lutions are the same, is not yet in the 6th-order-convergence 
regime, and for these choices of grid resolutions, appears to 
be 2nd-order convergent. We emphasize that this is not a 
demonstration that the simulations have entered a 2nd-order- 
convergence regime; while it is expected that at sufficiently 
high resolutions the 2nd-order-accurate components of the 
code will dominate, we have not yet performed any simula- 
tions of any configuration with high enough resolution to see 
clean asymptotic 2nd-order convergence. All the second panel 
of Fig. [Tltells us is that we are not yet in the fully convergent 
regime, but since the results are converging, in the sense that 
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FIG. 1. Convergence of the phase as a function of time for the %[ = 0.5 (left) and q = A (second convergence series) cases. The early-time 
behaviour is shown in the upper plots, and the late-time behaviour in the lower plots. Scaling with respect to different convergence orders 
is shown, to illustrate how cleanly the data exhibit a particular convergence behaviour. In these plots / = indicates the beginning of the 
simulation, and A(|)(? = 0) = in all simulations. The Xi = 0-5 case shows reasonably clean 6th-order convergence, and the accumulated phase 
difference is A0 = 0.43 rad between the medium- and high-resolution simulations. The g = 4 case is not yet in the 6th-order convergent regime, 
and appears (erroneously) to exhibit 2nd-order convergence. The accumulated phase difference between the medium- and high-resolution 
simulations is 1.5 rad. 



the eiTors reduce between simulations, we can still make a 
conservative estimate of the accumulated phase error 

In the Xi — 0.5 case, we can use Richardson extrapolation 
based on 6th-order convergence, to estimate the uncertainty 
in the accumulated phase as 0.6 rad. In the q — \ case, how- 
ever, we do not yet see 6th-order convergence. If we were to 
optimistically assume that the medium- and high-resolution 
simulations are in the convergent regime, and it is only the 
low-resolution simulation that is not, then we would estimate 
an accumulated phase uncertainty of 1 .9 rad based on Richard- 
son extrapolation. If we instead produce a much more conser- 
vative error estimate based on 2nd-order Richardson extrapo- 
lation, we find 6.7 rad. Note that this is an order of magnitude 
higher than we found in the cleaner Xi = 0-5 case, and coiTe- 
sponds to a full GW cycle. 

How are we to interpret these accumulated phase errors? 
They are certainly useful in comparing simulations — for ex- 
ample, the second q — 4- convergence series is less accurate 
than the Xi = 0.5 seiies. But this measure of the phase ac- 
curacy is of little additional value. In a GW application we 
will use only the waveform after the passage of the pulse of 
junk radiation. It is then difficult to estimate the phase error 
of the resulting waveform because to do so we must first align 
the waveform at some point after the beginning of the simu- 



lation. This introduces ambiguities (due to numerical noise in 
the GW phase and frequency) that may be larger than the error 
we ultimately want to measure; this point is illustrated well in 
Sec. V.E.2of |68|. 

When we compare with PN approximants in Sec. |VII[ we 
will be interested in the phase eiTor over the 8-10 GW cycles 
up to GW frequency Mco,„ =0.1. If we look at Fig. [T] we find 
that the accumulated phase difference at that frequency, be- 
tween the medium- and high-resolution simulations, is 0. 1 rad 
and 0.12 rad, respectively for the Xi — 0.5 and q — 4- cases. 
However, if we instead Une up the waveforms in each con- 
vergence series at Mo),,, =0.1, and measure the accumulated 
phase disagreement as we go back 10 cycles, we instead find 
about 0.01 rad for both configurations. This is an order of 
magnitude lower than what we observe when the waveforms 
are aligned at the beginning. This is an artifact of both the 
removal of the junk-radiation portion of the waveform, and 
simply the properties of the waveform frequency functions; 
similar effects are seen with different choices of alignment of 
PN waveforms, for which no junk radiation or significant nu- 
merical noise exist. 

These results demonstrate that we must be careful to choose 
our assessment of the phase error consistently with the ap- 
plication we are interested in. For the PN phase comparison 
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in Sec. VII we compare PN and NR waveforms aligned at 
Mco,n — 0.1, and so the only meaningful numerical phase er- 
ror estimate that makes sense is that based on the same form 
of phase alignment. 

Noise in the numerical frequency introduces an ambiguity 
into the matching time for any phase re-alignment procedure, 
which makes it impossible to use the realigned phase as the 
basis of a convergence test. However, we can vary the match- 
ing time within its error bounds, measure the maximum ac- 
cumulated phase disagreement that arises from this process, 
and then use 4th-order Richardson extrapolation to provide a 
conservative error estimate in the phase. The results of this 



process are shown in Tab. Ill and will be relevant to the anal- 



ysis in Sec. |VII| The same procedure and alignment are used 
to give estimates of the phase uncertainty accumulated dur- 
ing merger and ringdown. The table also shows an estimate 
of the total accumulated phase error, based on a convergence 
analysis like that shown in Fig.[T] we repeat that this estimate 
has no direct relation to any physical application, and is only 
useful as a means to compare the relative accuracy of different 
simulations. Note that this number is not simply the sum of 
the inspiral and merger phase uncertainty estimates, and this is 
a clear artifact of the alignment ambiguity in assessing phase 
accuracy. As such, in most physical applications, where some 
realignment is implicitly performed, the effective total phase 
error may drop by an order of magnitude over the numbers 
shown in the table. It is also clear that the total accumulated 
phase error estimates depend dramatically on whether we see 
clean convergence (the one truly clean case shown in the table 
is Xi ~ 0.5; other lower spin cases are also cleanly conver- 
gent). Nonetheless, we will see in Sec. V C that this level of 



accuracy is still well within the requirements for GW detec- 
tion. 



B. GW amplitude 

We now consider the GW amplitude. This plays a less im- 
portant role in (ietection, but errors in the amplitude (as well 
as higher harmonics) will affect estimates of the source pa- 
rameters, since all parameter errors scale with inverse signal- 
to-noise ratio (SNR). 

If we perform a time-domain convergence analysis of the 
GW amplitude, our conclusions are biased because the appar- 
ent amplitude error is in fact a combination of both the ampli- 
tude and phase errors — if the amplitude were measured with 
no error by the code, but two waveforms are out of phase, they 
will appear to have a non-zero amplitude error when com- 
pared in the time domain. We discussed this point in some 
detail in [14K and used a parametrization of the amplitude in 
terms of GW phase to reduce the effects of dephasing on the 
amplitude analysis. This works well if the phase error as a 
function of GW frequency is small, but this will not always be 
true. We expect (from the PN and perturbation theories) that 
the GW amplitude is a function of the GW frequency, and so 
the ideal method to measure the amplitude accuracy would be 
to reparametrize the amplitude as a function of GW frequency. 

This procedure also presents problems: the GW frequency 



is a numerically noisy function during the early and late parts 
of the simulation; it is certainly not the smooth monotoni- 
cally increasing function that we expect it to be on physical 
grounds. We can partially circumvent this difficulty by pro- 
ducing a smooth analytic fit of the frequency function, and 
considering the GW phase and amplitude as parametrized by 
that function. The smoothing process may itself introduce nu- 
merical artifacts, and either mask or exaggerate the conver- 
gence properties of the numerical results. But in general it is 
sufficient to allow us to calculate uncertainty estimates for our 
waveforms. 

Our method for modeling the GW frequency is as follows, 
based on an earlier version that was used (for equal-mass, non- 
spinning waveforms) in the work for the Samurai project 1 1 3 1 . 
For the inspiral, we start with the analytic TaylorT3 approx- 
imant for the frequency, as given in |69|. We neglect the 
highest-order (3.5PN) nonspinning term and replace it by a 
free parameter that will be fit to our data. In addition, follow- 
ing 1.69.1 , we do not specify the value of the spin, but also treat 
it as a free parameter — remember that our goal is to produce a 
clean analytic fit to the frequency, and we are not interested in 
whether all of the parameters have their usual physical inter- 
pretation. The modified TaylorT3 frequency function is then 
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where v — mim2/M^ is the symmetric mass ratio, S = Si +S2 
is the total spin parallel to the orbital angular momentum, 
E ~ M{S2/ni2 — S\ /mi), and 5M — mi— m2- (Note that in the 
Samurai paper, the PN frequency formula Eqn. 7 is missing an 
overall factor of two.) In the cases we consider here, the spins 
are nonzero only in the equal-mass case, and the spins are 
always equal to each other, so the (5ME) terms do not con- 
tribute. The function T is usually given by T = v{tc — t)/{5M), 
and tc is interpreted as the "time of coalescence" in standard 
PN theory, although a more appropriate term would be "time 
of divergence". 

In order to produce a formula that can be fit through our 



data, we redefine T as 



25M2 



(2) 



where both tc and d are free parameters that are fit to the data. 
This modification of T prevents HpN from diverging at f = tc- 
In the form that we have written it, i2pN is now symmetric 
about t — tc, which is certainly not physically realistic, but be- 
yond this point we will make a smooth transition to a different 
function, which models the ringdown. 

To model the ringdown phase, we modify the ansatz sug- 
gested in llTOl . and write the full frequency as 



n(f) = iipN(T)- 



[%-i2pN(T)] 



1 + tanh[ln ^/K-{t- k))/b\ 



(3) 



The constants {fc,fo,5, K^a^b^Hf} are parameters that are de- 
termined to produce the best fit to the numerical data. The 
constant £lf corresponds to a fit of the ringdown frequency, 
but the other parameters have no clear physical interpretation. 
(Even the "spin" parameter S really amounts to no more than a 
modification of the 1 .5PN and 2.5PN terms in the description 
of the inspiral frequency.) 

Fig.l2]shows a typical frequency fit, in this case for Xi = 0-5. 
We see that the dominant error in the fit is due to the resid- 
ual eccentricity in this simulation; recall that the aligned-spin 
cases are based on QC parameters and have the highest ec- 
centricity of all the cases we studied. The procedure does not 
work quite so well in cases with high spin; the frequency evo- 
lution is not captured so well during the early inspiral, or in 
the 2QQM before the peak GW amplitude. The fitting formula 
(|2]i could be modified to address this, and indeed the model of 
the transition to ringdown ^ has since been improved by the 
authors of ITTOIItTI . These issues, and the masking of eccen- 
tricity effects, mean that this frequency fit is far from ideal, 
and cannot be used for a convergence study of the amplitude. 
However, it is adequate for the purpose of providing a rough 
estimate of the amplitude uncertainty in our simulations. 

Fig. [3] shows the differences in A{(o) with respect to res- 
olution for the Xi — 0.5 configuration. The figure suggests 
that the error in the GW amplitude due to numerical resolu- 
tion is on the order of 1%. At late times the relative error 
grows higher, but this is beyond the frequency at which the 
amplitude reaches a maximum (indicated by the dashed line), 
and is well into the ringdown of the signal. Note that if we 
perform an error analysis based on the time-domain ampli- 
tude, then the maximum error between the medium- and high- 
resolution simulations is around 6%, which suggests that it is 
indeed dominated by the phase error 

Estimates for the amplitude uncertainty (using the ampli- 
tude parametrized by GW frequency) are given in Tab. Ill In 



all cases the GW signal was extracted at 7?ex = 90M; we will 
discuss the errors due to the use of a finite extraction radius in 
the next section. 



C. Mismatch with respect to numerical frequency and GW 
extraction radii 

Ultimately we are interested in the accuracy of our wave- 
forms with respect to GW detection. The most meaningful 
way to do this is to calculate the faithfulness between wave- 
forms from different numerical resolutions and different ex- 
traction radii. 

A calculation of the faithfulness is based on the overlap be- 
tween two waveforms. The overlap is usually calculated in 
the frequency domain. For two GW signals (in this analy- 
sis all quantities are with respect to the (i — 2,m — 2) mode 
of the signal) /ii(/) and hiif), we define an inner product 
weighted by the power spectral density of the detector noise, 
5„(/),as||72l, 
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Our data represent ^4(f), not the wave strain h{t), but the 
two are related by ^4 — h+ — ihx- Making two time integra- 
tions is trivial when transforming to the frequency domain, 
and although this does not automatically remove the irrita- 
tion of having to choose constants of integration (see 1.73] for 
a recent discussion of this problem) we have found that our 
ignorance of these constants does not affect mismatch calcu- 
lations 113|. 

Given the definition of the inner product {hi\h2), we nor- 
malize it and maximize over phase and time offsets in the 
data. If the waveforms were equal, then this quantity would be 
unity. This is the faithfulness of the waveform: it is a measure 
of how "far" a theoretical waveform is from a supposedly true 
waveform with the same physical parameters. We define the 
faithfulness mismatch as the deviation from unity: 



1 — max — , , 

T,* V(/2l|/2l)(/22|/22) 



(5) 



Ideally the integrations over frequency are in the range 
[0,°°]. When we have a finite data set like our numerical wave- 
forms, we also need to optimize with respect to the window 
of our data that we sample. This is discussed in farther de- 
tail in 1 13 1; the optimization with respect to phase and time 
offset is trivial in the frequency domain when using only one 
the {i = 2,m = 2) mode. In general more sophisticated max- 
imisation procedures are required, for example the techniques 
described in |74| for the quadmpole harmonic, and in llTSJI for 
signals that include higher harmonics. 

The faithfulness mismatch is calculated without any opti- 
mization over the intrinsic parameters of the binary. In a true 
GW search using a bank of theoretical templates, one opti- 
mizes not only over time and phase shifts, but over all phys- 
ical parameters included in the template bank. Optimization 
over the physical (intrinsic) parameters gives the effectualness 
mismatch, i.e., how well a waveform family will be able to de- 
tect GW signals irrespective of whether the physical parame- 
ters are measured correctly. We have access to waveforms 
representing only one choice of intrinsic parameters of the bi- 
nary, and so cannot perform this optimization (although we 
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FIG. 2. Analytic fit to the GW frequency for the Xi ~ 0.5 case. The right panel shows the fractional difference between the fit and numerical 
data. For this configuration, the error in the fit is dominated by the residual eccentricity in the simulation. The dashed line indicates the point 
at which the amplitude reaches its maximum. 

TABLE ni. Estimates of uncertainty in phase and amplitude. The phase uncertainty accumulated during the inspiral is based on an alignment 
of the GW phase at Ma = 0.1, and includes only the ten GW cycles up to that frequency, for consistency with the analysis in Sec. |VII| The 
same alignment is used for the phase uncertainty of the merger and ringdown regime. The "complete" phase uncertainty is a conservative 
estimate of the total accumulated phase error over the entire waveform, and is only relevant for relative comparisons of different simulations; 
see text in Sec.|V A| The amplitude uncertainties are described in Sec.|VB[ and the mismatch errors in Sec.|VC[ 
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could easily optimize over the total mass of the binary). The 
physical-parameter-optimized effectualness mismatch will al- 
ways be better than (or equal to) the faithfulness mismatch 
and so we can use the faithfulness to set an upper bound on 
the error of the waveforms. 

The faithfulness mismatch between the simulations of dif- 
ferent resolutions is negligible — it is below 10^^ for all rel- 
evant masses (down to about IOOMq) with respect to the Ad- 
vanced LIGO noise curve 1.76.1 (we use the approximate an- 
alytical formula displayed in ll36l ). where we choose a low- 
frequency cut-off of 20 Hz. 

We also use the faithfulness to estimate the error due to the 
finite extraction radii. The GW signal is extracted on spheres 
of radii /?ex = {50, 60, 70, 80, 90 }M, and we expect the error 
relative to the true signal as /?ex — > °° to fall off as l/R^x- This 
error is typically larger than that due to finite-difference er- 
rors, i.e., the finite extraction radius is the dominant source of 
error in the simulation. We estimate the mismatch error by 
extrapolating to R^^ -^ °o the mismatches between our finite- 
extraction-radii data. We find that the maximum mismatch 



(which is always at the lowest mass we consider, lOOM©) is 
2.8 X 10^'*. This is much larger than the mismatch due to the 
numerical resolution errors, as we expect. The maximum mis- 
match in each case is given in Tab. 



Ill 



Such levels of accuracy are well within the requirements 
set out in |77|, and also within the suggested accuracy for 
waveform modeling within the NR-AR project 1781 . This 
is also comparable or better than the level of mismatch be- 
tween the different equal-mass nonspinning waveforms (taken 
from independent codes) that were studied in the Samurai 
project 113], suggesting that these waveforms are also of 
sufficient accuracy for GW detection purposes with current 
ground-based detectors. The accuracy requirements for pa- 
rameter estimation, and for applications with future detec- 
tors, such as the space-based LISA |79| and third-generation 
ground-based Einstein Telescope |80|, may be much higher, 
but have not yet been quantified for NR waveforms. 
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FIG. 3. The amplitude error as a function of GW frequency for the 
Xi = 0.5 case. The deficiencies of the frequency fitting procedure 
preclude the use of A(co) for a convergence test, and the differences 
between the low-, medium- and high-resolution simulations are not 
scaled in any way. 



VI. PHYSICAL PROPERTIES OF THE BINARY 
CONFIGURATIONS 

Now that we have established the accuracy of our simula- 
tions, we can calculate some of their physical properties. The 
most accurate are quantities calculated from the phase and am- 
plitude of the leading harmonic, like the mass and spin of the 
final black hole. Less accurate are integrated quantities based 
on the leading subdominant harmonics, like the radiated en- 
ergy in each mode. The gravitational recoil, which is not only 
based on the higher harmonics, but on overlaps between some 
of the weaker harmonics, is the least accurate. We will con- 
sider first the general physical properties of the binary config- 
urations. 



A. General properties 

In Table|lV|we indicate the initial coordinate separation of 
the binary D/M, and the number of GW cycles before merger, 
Afew- The latter quantity is defined as A<I>/(27i;), where A<I> 
is the accumulated GW phase from t — 2QQM (i.e., after the 
early burst of junk radiation) until the time when the wave's 
amplitude reaches its maximum value. 

The configurations we simulated for this work clearly 
demonstrate the orbital hang-up and "anti-hangup" effects for 
spins parallel or anti-parallel to the orbital angular momen- 
tum. The orbital hang-up case was first studied in lISTl . and 
for a larger range of cases in |T6l; the largest spin considered 
was Xi ~ 0-92 in 1601 . One case of anti-parallel spins with 
Xi = -0.438 was considered in ||82| . 

When the black holes are nonspinning, a binary with an 
initial coordinate separation of D = 12M produces around 19 
GW cycles before merger When the black holes have spins 
X, = 0.25, the merger is delayed in comparison with the non- 
spinning case, and a binary with the same initial separation 
produces 21.5 GW cycles before merger Conversely, when 
the spins are Xi = —0.25, the merger is accelerated, and the 



binary produces only 18.5 GW cycles before merger. These 
trends continue as the spins are increased, and in order to 
produce comparable numbers of GW cycles in each simula- 
tion, the initial coordinate separation is increased for increas- 
ing anti-parallel spins Xi < 0, and decreased for increasing 
parallel spins Xi > 0. For the highest-spin cases, \Xi\ = 0.85, 
an initial separation of D = 13M is required to produce 16 
GW cycles in the anti-parallel case, and an initial separation 
of only D = lOM produces 20 GW cycles in the parallel case. 
In the unequal-mass nonspinning cases, we see that the 
number of cycles before merger also varies with the mass ra- 
tio q. The general effect is best understood by considering 
the two extreme cases, q— I and the extreme-mass-ratio case 
q^oo. In the extreme-mass-ratio case, i.e., a point particle or- 
biting a Schwarzschild black hole, the slow inspiral terminates 
abruptly at the innermost stable circular orbit (ISCO), and the 
small black hole plunges into the large black hole. Prior to the 
ISCO, the small black hole follows a slow adiabatic inspiral, 
so that there are very many orbits at separations just above the 
ISCO, but as soon as the small black hole passes the ISCO, 
there are no more orbits, only the fast plunge. By contrast, in 
the equal-mass nonspinning case, the transition from "inspi- 
ral" to "plunge" is very smooth, and there is no ISCO; the rate 
of inspiral simply increases. As the mass ratio is increased, the 
rate of the "plunge" increases, and the rate of inspiral prior to 
merger decreases — in other words, the dynamics approach 
the extreme-mass-ratio situation, and the system gets closer 
to exhibiting an ISCO. This behaviour is illustrated in Fig. p] 



B. Final mass and spin 

The final mass of the merged black hole can be estimated 
from the energy lost through gravitational radiation. Given the 
total (ADM) energy in the initial data, /^adm, and the radiated 
energy Eraj, we know that the final spacetime must contain the 
energy Mf = /sadm — £^rad- Since the final spacetime contains 
only a single stationary (i.e., non-radiating) Kerr black hole, 
Mf must be the mass of that black hole. 

We calculate the radiated energy E^a on each of the five ex- 
traction spheres R^^ = {50, 60, 70, 80, 90 }M, and extrapolate 
the result to /?ex — > °° assuming that the error falls of as 1/^ex- 
This assumption is most consistent with the data at the largest 
extraction radii, and so we use only 7?ex = {70,80,90}M for 
the fit, and include also R^^ = 60M to assess the robustness 
of the result. Once E,-^^ has been estimated for each resolu- 
tion, we find that the results converge at roughly 4th-order, 
although since the convergence is not extremely clean, we use 
2nd-, 4th- and 6th-order Richardson-extrapolated values to es- 
timate the uncertainty in the value from the highest-resolution 
simulation. In all cases we consider the uncertainty in the ra- 
diated energy to be about 2%. The values of the mass of the 



final black hole are given in Tab. IV 



To estimate the spin of the final black hole, we make use of 
analytic results that give the quasi-normal ringdown frequency 
MfCORD in terms of the black-hole spin, Of/Mf |83|. In the 
ringdown regime, the GW signal behaves as ^ exp(— /Ord?), 
where Ord consists of a real part (which is the frequency of 
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FIG. 4. Puncture motion for the nonspinning binary configurations q = I (left) and q = 4 (right). The figure shows about seven orbits before 
merger for each system. In the q = 4 case, the black line indicates the small black hole, while the red line indicates the large black hole. Note 
that the transition from inspiral to plunge and merger is more gradual in the q = I case. As the mass ratio increases, the plunge begins to 
resemble the ISCO effect that is present for extreme mass ratios. 



the ringdown waveform), and an imaginary part, which de- 
scribes the rate of exponential fall-off. Given the final mass 
Mf and the ringdown waveform, we can estimate the final 
spin using either the exponential decay-rate of the wave's am- 
plitude, or the wave's frequency M(0 in the ringdown stage. 
We find that matching to the ringdown frequency gives the 
most accurate results, in the sense that both methods agree 
within uncertainties, but the uncertainty estimates are smaller 
when we match to the ringdown frequency. In general our fi- 
nal spin estimates have an uncertainty of 1%, although it is a 
little smaller in the Xi =0.85 case, which we now consider in 
more detail. 

In the Xi = 0.85 case, we find that the ringdown frequency is 
McOrd = 0.769 ±0.001; see Fig.|5] The final mass is Mf/M = 
0.895 ±0.015, and the final spinis aj/Mj = 0.915 ±0.007. 
Note that the final mass is lower than quoted in llT6l . where all 
of the analysis was performed on the highest-resolution wave- 
form calculated on the largest radiation extraction sphere. 
Here we extrapolate the results with respect to extraction ra- 
dius (assuming a I /Rex fall-off in the error), and with respect 
to numerical resolution, where the results show between 2nd- 
and 6th-order convergence. The radiated energy increases 
with extraction radius, and so our estimate of the final mass 
decreases; this is why our extrapolated value (0.895) is lower 
than the /?ex = 90M value of 0. 9 1 1 quoted in 1 1 6 1 . In addition, 
we estimate the ringdown frequency using a 50M-long sample 
of the waveform starting 50M after the peak amplitude of the 
{£ = 2,m = 2) mode, while our earlier results were based on a 
portion of the waveform starting only a few M after the peak 
amplitude, which distorts the final estimate of the ringdown 
frequency. 

For comparison, Dain, et. al. Il60l study the Xi — 092 case. 



The initial black-hole spins are larger than studied here (and 
were set up to approach the highest spin possible for Bowen- 
York data), and therefore the final black holes should have a 
larger ringdown frequency than in our Xi = 0.85 case, and a 
higher final spin, and this is indeed the case. Note that they 
use the ADM mass Madm as the defining length scale in their 
simulations, while we use the total BH mass. In terms of the 
ADM mass, the ringdown frequency for the Xi = 0.85 case is 
A/ADMfita) = 0.761, while the frequency found in |60| for the 
Xi = 0.92 case is MadmOrd = 0.766. In addition, they give 
a measure of the final spin between 0.910 and 0.916, where 
0.915 is the value obtained using the same quasi-normal-mode 
method that we have applied here. Their final-spin result is 
consistent with ours' within our error bounds. 

A number of aligned-spin cases were studied in ['84 1. For 
the cases where direct comparison is available, our results 
show agreement within 1% for spins up to 0.5, and within 2% 
for higher spins. Assuming that their uncertainties are com- 
parable to ours', then the results agree. The first hang-up and 
anti -hang-up cases were studied in lISTl . each with spin values 
of \Xi\ = 0.757. They estimate final spins of 0.443 and 0.890 
for the anti-hangup and hangup cases respectively, and these 
are also consistent with our results. Two of the unequal-mass 
cases were also studied in fTOl, q = 2,4, and the final mass 
and spin results are in excellent agreement. 

We have also compared our results with fits for the final 
spin available in the literature. We find excellent agreement to 
about 1 % or better with f85l, as well as with |86| as long as 
the spins are not anti-aligned with the orbital angular momen- 
tum. For the latter paper we find disagreements of sa 10% for 
the cases Xi = -0.75, -0.85. 
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TABLE IV. Summary of the configurations simulated. The table indicates the initial coordinate separation D/M of the punctures, their 
tangential and radial momenta {pi,p,-), and the eccentricity e of the resulting coordinate motion. For the Xi < cases, where enhanced PN 
parameters were used to achieve low-eccentricity inspiral, the eccentricity from raw PN-inspiral parameters is also shown in brackets. The 
initial GW frequency is Mft),, and the ringdown frequency of the final merged black hole is Mcord- The simulation includes Wgw cycles before 
the peak of the GW amplitude, which occurs at fpeak- The final black hole has mass My and spin ay, and receives a recoil of vi^jck- 
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Pt/M 


-p,/M(xlO-4) 
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A'gw 
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Gf/Mf 


Vkick (km/s) 
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0.0025 (0.009) 


0.040 


16 


1868 


0.457 


0.969 


0.412 
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0.049 


20 


1739 


0.650 
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FIG. 5. The numerical GW frequency at a time {t — fpeak)/^ af- 
ter the peak of |r1'4 22|! shown for simulations at three resolutions. 
The frequency oscillates around a value that we take to be the ring- 
down frequency. The amplitude of the oscillations decreases as the 
numerical resolution is improved, suggesting that these are only a nu- 
merical artifact. The upper and lower bounds of our estimate of the 
ringdown frequency are indicated by the two horizontal lines in the 
plot. The bounds were obtained by considering the results from all 
three numerical resolutions, and varying the portion of the data used 
for the fit; the final quoted values were calculated using the range 
{t — fpeak)/^ £ {50, 100}. Despite the high amplitude of the noise 
in the data, the average value shows very little variation, and can be 
estimated with an uncertainty of only Aco = 0.001/M. 



ifTOl El). However, knowledge of the subdominant modes 
may aid detection, and are important for accurate estimation 
of the source parameters ll64l[89]492ll . 

We assess the relative importance of the sub-dominant 
modes by calculating the energy radiated in each mode. The 
radiated energy in each mode is given by | 



—-— = hm -—- 



^ilmdt 



(6) 



In practice the limits of the integration are taken as the time 
in the simulation just after the junk radiation has passed, and 
a time after the signal has rung down to the level of numerical 
noise. The results are summarized in Tab. IV] including only 
those modes that contribute above 1 % of the total energy. We 
see that in the equal-mass cases, the (2,2) mode dominates — 
around 98% of the energy is radiated in the dominant mode in 
all cases, with only a negligible variation due to spin, and no 
other modes contribute above 1%. 

In the unequal-mass cases, the energy contribution from the 
higher harmonics grows rapidly with mass ratio, particularly 
in the £ — zLm modes. We defer the reader to the detailed 
discussion in |94|, but note that, even at q = 4, most of the 
energy is radiated in a very small number of harmonics. 



D. Recoil 



C. Energy spectrum In spherical harmonic modes 

The {£ — 2,m = 2) mode dominates the GW signal from 
a black-hole-binary coalescence, and indeed most current 
searches in detector data employ templates that include only 
this harmonic |87 88| (see also the NINJA project searches in 
simulated data with injected numerical relativity waveforms 



Due to the asymmetry of the radiation emission in the 
unequal-mass cases, linear momentum is radiated from the 
system, and the center-of-mass of the binary moves as the 
black holes inspiral. The direction of the center-of-mass recoil 
rotates with the binary, so that the average movement is small. 
However, the rate of momentum loss grows as the black holes 
get closer, and, as with the total GW signal, peaks at merger 
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TABLE V. Ratio of total energy radiated in each mode. Only contributions above 1% are included. 



Case 


(2,±2) 


(2,±1) 


(3,±3) 


(4, ±4) 


(5,±5) 


Xi = -0.85 


0.988 








— 





Xi = -0.50 


0.989 








— 





Xi-0 


0.990 








— 





Xi^Q.5Q 


0.988 








— 





;t; = 0.85 


0.988 








— 





q = 2 


0.947 


— 


0.038 


— 


— 


q = 3 


0.897 


— 


0.076 


0.013 


— 


q = A 


0.868 


0.013 


0.095 


0.017 


— 



This final burst of GW emission causes an overall recoil, or 
"kick". 

Since the bulk of the recoil arises during the merger, short 
simulations are sufficient to accurately measure the effect, 
and these were used in |25 95, 96 1 to make the first accu- 
rate fully general-relativistic predictions of gravitational re- 
coil, and found that the maximum kick for nonspinning bina- 
ries is Vniax = 175 ±11 km/s for a mass ratio of q = 2.8 ll25l . 
An analytical fitting formula for the recoil from nonspinning 
binaries was presented in |25|, for recent papers containing 
such fitting formulas see |97| (which uses the same ansatz as 
1,25 1 and finds slightly different but consistent fitting parame- 
ters) and L86J (which quotes the fit from 1251 for nonspinning 
binaries). 

Tab. |IV] shows the results for the current simulations, 
which agree with those from the shorter simulations pre- 
sented in [25 1. It has also been shown that much larger 
recoils are possible from spinning or highly elliptical bina- 
ries |26, 60 98-1051, but not for any of the configurations 
that we have studied in this work. 

Our newer simulations improve over those produced in |l25l 
in two ways: they include many more cycles before merger, 
and the wave extraction is performed at larger radii. On the 
other hand, the numerical resolution at the wave extraction 
radii is lower, which reduces the accuracy. As such, the values 



we quote in Tab. IV have large error bars. 

The flux of angular momentum radiation is given by 



dPi .. 

—— = lim 

at r^'" 



r 
16n Jn 



^'4dt 



dQ. 



(7) 



where £{ ~ (sin0cos0,sin0sin^,cos0) 01O6L The total re- 
coil is calculated by integrating Eqn. (|7]) over the duration of 
the simulation. 

The additional length of the new simulations allows us to 
remove one source of error in shorter simulations: the choice 
of starting time in the integration of dPi/dt to calculate the 
total radiated linear momentum. This function oscillates with 
time, and during the inspiral the average radiated linear mo- 
mentum is much smaller than the amplitude of the oscillations 
— so a poor choice of starting time in the integration of dPi/dt 
could potentially corrupt the final result. In |25| the uncer- 
tainty due to this effect was estimated at about 3%. In |[95ll98l 
attempts were made to both account for this effect and for the 



linear momentum loss that will have accumulated over the ear- 
lier inspiral of the binary. In our cases, where we possess the 
waveform for many more cycles before merger, we are able 
to simply calculate the total recoil for a range of integration 
starting times to, and then to take the average of these values. 
We find that the uncertainty in this process is only a fraction 
of a percent of the final result. 

Figurel6]illustrates this effect with the q = 4- case. The lower 
limit of the integration, Iq, was varied between Iq — 125M (just 
after the burst of junk radiation has passed through the signal, 
and to = 1200M, which is about 300M before dPi/dt has fallen 
to negligible values, and also roughly corresponds to the value 
of fo that was used for the much shorter waveforms studied 
in |25 1. Note that the kick calculated for different choices of 
to varies by about 4 km/s, or 3% of the result. A linear curve fit 
through the results (shown in the figure) indicates that the av- 
erage result of the integrated linear momentum radiation rises 
very slowly during the inspiral, and varies by only 0.7 km/s, 
or 0.5%. In our results, we determine the final kick to be the 
average over this range of choices fo, which introduces only a 
negligible error in our result. 

A second error source that could not be quantified in ||25]| 
was that due to extraction of the GW signal at finite extraction 
radius. In that work an extraction radius of only /?ex = 30M 
was feasible. We now extract GW signals at up to Rex — 9QM, 
although we find that the numerical resolution at the wave ex- 
traction spheres allows an accurate calculation of the recoil 
only for R^x = {50,60}M. However, these two radii are suf- 
ficient for us to extrapolate the recoil to /?ex -^ °°, assuming 
a 1 /Rex fall-off in the error. This gives the values listed in 



Tab. IV This fit also of course predicts the value of the recoil 
at Rex = 30M, which agrees well with the values in 1 25 ] . How- 
ever, due to the poorer numerical resolution on the extraction 
spheres, we assign large error bars to our values. 



VII. LATE INSPIRAL COMPARISON OF NR AND PN 
WAVEFORMS 

One of the most important applications of our waveforms 
is as input in the construction of analytic waveform mod- 
els that can in turn be used to construct template banks for 
GW searches. In particular, these waveforms have already 
been used to produce the phenomenological models presented 
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FIG. 6. Variation in the estimate of the total radiated linear momen- 
tum (recoil), as a function of the starting time fo of the integration of 
dPi/dt for the q = A case. The value oscillates around a slowly grow- 
ing average, which is indicated by a straight line. Note that the latest 
integration time used, fo = 1200M, is approximately 300M before 
the end of the ringdown, which is the point at which the integration 
was started in the older calculations in 1251. 



in ll36l l56l l65l l67l . in which the NR waveforms for the late 
inspiral and merger are connected to long PN inspiral wave- 
forms, to produce "complete" waveforms for the full inspiral- 
merger-ringdown, and it is these complete waveforms that are 
then used in the construction of a phenomenological model 
that is essentially an analytic fit across the relevant section of 
the black-hole-binary parameter space. 

In performing this procedure, we need to quantify the level 
of agreement between the PN and NR waveforms in some 
region where they are both considered to be valid. In other 
words: the approximate PN waveforms are expected to be ac- 
curate during much of the long inspiral, but are they still ac- 
curate enough at the point where we want to connect them to 
fully general relativistic results? 

This question was first addressed for equal-mass nonspin- 
ning waveforms in 169^, '1071, and later with increasing levels 
of precision in |14, 15, 68 1. The conclusion of these works 
was that the phase disagreement between NR waveforms and 
typical PN approximants was less than 1 rad over the last 10 
cycles up to Ma = 0.1, and the error in the quadrupole PN 
amplitude was about 8%; however, the phase could be tracked 
with surprising accuracy by one PN approximant, TayIorT4, 
and the PN amplitude was accurate to within 2% when evalu- 
ated at 3PN order 1681 [T08l. 

Similar comparisons were performed with the equal-mass 
Xi > cases that we consider here, where it was found that 
the phase disagreements were comparable for all spin values 
for the TaylorTl approximant (i.e., approximately 1 rad over 
the 10 cycles up to Mco = 0.1), but that the TayIorT4 approx- 
imant, which performed so well in the nonspinning case, did 
no better than TaylorTl (and often worse) when the BHs were 
spinning. It should be noted, however, that the Taylor approx- 
imants did not include spin terms up to the same PN order 
as in nonspinning terms (2.5PN versus 3.5PN), which is a 
point we will return to later In addition, it was found that the 
quadrupole amplitude error grew to as much as 12% in high 



spin cases fl6]. Equal-mass nonspinning eccentric binaries 
were considered in |109|, and one unequal-mass precessing- 
spin configuration was studied in ||37]| . 

In this section we will perform a PN-NR comparison for all 
of our waveforms, which now include ^ ^ 1 and Xi < cases. 



A. PN approximants 

The PN approximants considered here are derived from the 
energy S" and GW flux ^ of a black-hole binary on quasi- 
circular orbits. Both quantities are given in the PN framework 
as expansions in v/c, up to [v/cf (3.5PN order), where v is 
the relative velocity and c the speed of light. Following the 
standard convention, we regard S" and ^ as functions of the 
dimensionless variable x = (v/c)^ that is related to the orbital 
phase 0orb via 



X = (MOJorb) 



2/3 



d0, 



'orb 



dt 



= Oorb 



(8) 



The energy-balance law dS" /dt = — ^ can be transformed to 
an evolution equation for x. 



dx 

dt 



^ 



dS' /dx 



(9) 



which in turn leads to the bn-moAc of the gravitational wave 
strain 



t'm ( 



h'"\t)^H'"\t)e 



-'m<l>oibi') 



(10) 



The amplitudes //^'" are given as expansions in x to 3PN order 
in the non-spinning case 1 1 10 1 and up to 2PN order in spinning 
contributions |59|. 

A direct (numerical) integration of (|9]) and ^ is referred 
to as the TaylorTl approximant. If instead the right-hand side 
of Eq. (|9]l is re-expanded as a Taylor series in x before inte- 
grating, the resulting approximant is called TaylorT4. This 
re-expansion is truncated at the same order as the energy and 
flux (i.e., 3.5PN); all higher powers in x are incomplete and 
therefore neglected. 

If we apply the same strategy to the spin contributions that 
enter at 1 .5PN (leading order spin-orbit coupling), 2PN (spin- 
spin) and 2.5PN order (next to leading order spin-orbit), we 
should neglect all spin-dependent terms in the re-expansion of 
Q that appear at 3PN and 3.5PN order We denote the result- 
ing approximant that was used for instance in [ 16] as TaylorT4 
(truncated). If we instead disregard the distinction of spinning 
and non-spinning terms and use the "full" re-expansion up to 
3.5 PN order, thereby keeping incomplete spin contributions 
at 3 and 3.5PN order, we denote the resulting approximant 
simply as TaylorT4. For a detailed discussion and explicit ex- 
pressions for the approximants see 1671 and references therein. 
Further small corrections to the spin contributions to the PN 
phase and amplitude, due to typographical or other eiTors in 
the original literature, were found during a program of PN- 
approximant verification within the Ninja collaboration [78|; 
these will be described in more detail in an upcoming amend- 
ment to the data format specification document lllllL and are 
discussed further in Sec. lVIICl below. 
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B. Phase comparison 

We now compare the PN and NR phase. Our procedure, as 
in previous studies [ 14, 16 1, is to consider the phase for the A^ 
GW cycles up to the matching frequency Maim =0.1. We line 
up the PN and NR phase functions so that they agree when 
(B = 0),,,, and relabel this event as f = 0. We then calculate 
the phase disagreement as it accumulates over A^ cycles back 
in time. Note that although our comparison is over a fixed 
number of GW cycles, it is not over a fixed frequency range, 
due to the different frequency evolution in each configuration. 
In the same way, the comparison is also over different lengths 
of time between different configurations. However, we have 
found that the qualitative behaviour of the comparison results 
does not depend on whether we compare over a fixed range of 
cycles, frequency, or time. 

In previous studies we simply calculated the phase differ- 
ence A0(f) = 0pn(O^^r(O' and quoted A0(fAr) as the accu- 
mulated phase difference, where t^ is the time A^ cycles prior 
to the point where (O — O,,,. This procedure gives consistent 
results, but we may worry in general that A0(r) is not a mono- 
tonic function, and so a more robust procedure is to consider 
instead 



A0(fw) = 



1 



0. .2 ll/2 

0NR(f)-0PN(O) dt 



tN 



(11) 



This gives us a measure of the average rate of increase of 
the phase disagreement. A similar procedure was also used 
although in that study the alignment of the wave- 



in IITT2I 

forms was adjusted to minimize A0. An elegant alternative 
measure of the accumulated phase disagreement is given in 
Eqn. (3.15) of 111131 . We instead wish to evaluate how well 
the PN phase evolution agrees with the fully general relativis- 
tic NR results. For comparison with previous results in the 
literature, we will also show the results of a direct calculation 

of 0pN(O-^R(O■ 
Fig. [t] shows the disagreement between the PN and NR 
phase for the equal-mass configurations with non-precessing 
spins over A^ = 10 GW cycles. Three PN approximants are 
used: TaylorTl, TaylorT4, and TaylorT4-truncated, as de- 
scribed in the previous section. 

We see that in both calculations of the accumulated phase 
disagreement, TaylorTl is the most robust. It performs best in 
the nonspinning case (which is to be expected, since the non- 
spinning contributions are known to higher PN order than the 
spinning contributions), and for all spinning cases the accu- 
mulated phase disagreement is between 1.0 and 2.0 rad, while 
the square-averaged phase disagreement is between 0.5 and 
1.0 rad. We see also that TaylorT4-truncated performs worse 
as the spin is increased, and for large anti-aligned spins per- 
forms very poorly. The full TaylorT4 approximant performs 
better for most spin values, although it is again poor for large 
anti-aligned spins. It is in light of comparisons using only 
TaylorTl and TaylorT4-truncated that we chose to use the 
TaylorTl approximant in the construction of hybrid wave- 
forms for the phenomenological model in ll56l . 

Fig. [8] shows a similar plot, but this time for the unequal- 
mass nonspinning configurations. The q — 2 simulations con- 
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FIG. 7. Phase disagreement between NR and PN results for three 
choices of PN approximant, for configurations that consist of equal- 
mass binaries with equal spins oriented parallel or anti-parallel to the 
orbital angular momentum. The first panel shows the accumulated 
phase disagreement for the ten GW cycles up to Mo>i„ =0.1. The 
second panel shows the integrated square of the phase disagreement, 
Eqn. flT). 



sist of less than ten cycles before Mo = 0.1, so we consider 
only A^ = 8 cycles in the phase comparison. In this case we 
see that TaylorT4 continues to perform well for unequal-mass 
configurations. We expect that at higher mass ratios the per- 
formance of all PN approximants will deteriorate, but up to 
(7 = 4 this deterioration cannot be clearly measured; the per- 
formance of TaylorTl and TaylorT4 shows some variation 
with mass ratio, but this is not monotonic. 

From our phase comparison analysis, we conclude that the 
TaylorTl approximant is most robust over the entire subset 
of the black-hole-binary parameter space that we have stud- 
ied. The TaylorT4 approximant performs well for all non- 
spinning cases. The performance of TaylorT4 for spinning 
cases varies greatly between our two choices of treatment of 
the higher-order spin contributions, but for both choices shows 
poor agreement for large anti-aligned spins. We caution, how- 
ever, that the performance of the approximants over a rela- 
tively small number of numerical cycles does not tell us how 
well they perform before at lower frequencies, and we will 
return to this point in the Discussion. 
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FIG. 8. Phase disagreement between NR and PN results for two 
choices of PN approximant, for configurations that consist of non- 
spinning black holes of unequal mass, with mass ratio q = M2/M1 . 



C. Amplitude comparison 

We now compare the PN prediction for the inspiral wave 
amplitude with numerical results, for the (£ = 2, m = 2) mode. 
We found in 1 14] that in the equal-mass nonspinning case the 
quadrupole PN amplitude was larger than the full GR ampli- 
tude during inspiral by about 7%. It was later shown in |681 
that the amplitude agreement could be improved to within 2% 
if corrections up to 3PN order were used. For equal-mass bi- 
naries with aligned spins, we found in fT6l that the quadiTipole 
PN amplitude disagreement rose to about 12% in highly spin- 
ning cases. 

In this section we extend our previous analysis of the 
quadrupole amplitude to anti-aligned and unequal-mass cases. 
We also compare with the PN amplitude that results from us- 
ing all cuiTently known amplitude corrections (up to 3PN or- 
der non-spinning M108I IllOl and up to 2PN order spinning 
contributions [59", 114]). We have taken care when combin- 
ing results for amplitude functions from different sources in 
the literature, in particular regarding different conventions for 
the choice of relative phase factors. In our implementation 
we now follow the convention of |59 1, which differs from that 
of miOL from which we originally took our nonspinning am- 
plitude contributions. We have checked for consistency with 
the amplitude of the I = \m\ = 2 modes as given in 11151 . and 
we have compared with an independent code as part of the 
Ninja project 178111161 . In addition, we have also checked 



that inconsistent choices of the relative phase factors (e.g., 
caused by misprints in the literature) significantly increase the 
deviation of the NR and PN amplitudes; the correct choices 
lead to the best agreement with results from full general rela- 
tivity. 

We find that the GW amplitude shows variations with nu- 
merical extraction radius that are comparable to the level of 
disagreement with the PN predictions. However, the error in 
the amplitude seems to fall off as l/R^-^ (see fT4\ for a dis- 
cussion of this effect), and allows us to perform an accurate 
extrapolation to /?ex —^ °°- Having obtained the accurate am- 
plitude of Rgy^'4, we then express the amplitude as a function 



of frequency, using the methods we introduced in Sec. V B 



which then allows us to easily compare with the PN ampli- 
tude, which is always expressed as a function of frequency. 
Note that for this comparison we perform a frequency fit to our 
data during only the inspiral, which allows us to much more 
accurately capture the amplitude evolution; it is now much 



more necessary than in Sec. V B to have a reliable physical 
fit. 

Fig. I9] shows the average disagreement between the PN 
and NR amplitudes over the 10 cycles up to Mco = 0.1, for 
the equal-mass spinning cases. The results using both the 
quadrupole and 3PN order amplitudes are shown. As seen 
in [161 the quadrupole amplitude disagreement rises to just 
over 12% for the highly spinning cases. The increase in dis- 
agreement is approximately linear with respect to the spin, 
and we predict that the maximum disagreement for extreme- 
spin black holes would be around 14%. For large anti-parallel 
spins, the quadrupole amplitude performs much better, and 
drops to around 3% for ;i^,- = — 0.85. 

When PN amplitude contributions up to 3PN (non- 
spinning)/2PN (spinning) order are used, the agreement with 
NR results is much better. In the nonspinning case it is 3%, 
consistent with the results in pMI. (Note that the uncertainty 
in the extrapolated NR amplitude is around 1 %.) The variation 
with spin is small, rising to only 4% in the high-spin hang-up 
cases, and falling to 2.5% in the high-spin anti-hang-up cases. 
We find similar results for the unequal-mass cases, where the 
average disagreement is around 3%. 



DISCUSSION 

We have presented the results of two sets of numerical 
simulations of black-hole binaries, equal-mass binaries with 
equal, non-precessing spins with Xi — St/Mf g [—0.85,0.85], 
and nonspinning unequal-mass binaries with q = M2/M1 G 
[1,4]. These simulations cover between six and 10 orbits be- 
fore merger. The most accurate simulations have a numerical 
phase uncertainty during inspiral of 0.05 rad, and a total accu- 
mulated phase error of about 1 .0 rad. The phase uncertainties 
in the least accurate case are 0. 1 rad during inspiral, and a to- 
tal accumulated phase error of up to 15 rad. We have shown, 
however, that the uncertainty estimates depend strongly on the 
alignment of the waveforms, and whether the results are repre- 
sented as functions of time or of GW frequency. The accuracy 
of the amplitude of the {£ = 2,m = 2) mode of ^4 is in gen- 
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FIG. 9. Average amplitude disagreement between PN and NR re- 
sults, over the last ten cycles up to Mco = 0.1. The quadrupole PN 
amplitude error is only about 3% for large anti-aligned spins, but 
rises to around 13% for large aligned spins. When the amplitude cor- 
rections are included up to 3PN-order, the PN amplitude error is only 
3-4% for all spin values. 



eral better than 1% during inspiral, and between 2% and 5% 
during merger. 

For purposes of GW detection, the important quantity to 
consider is the mismatch-error in the waveform. This is dom- 
inated by the errors due to the extraction of the GW signal 
at a finite radius from the source. However, in all cases the 
mismatch error (minimized over only time and phase) is be- 
low 10^^, meaning that the numerical waveforms are well 
within the accuracy requirements for detection with current 
and planned ground-based detectors. 

These statements of waveform accuracy for detection apply 
only to the dominant mode and, more importantly, are only 
relevant when we consider binary masses such that the en- 
tire numerical waveform is within the sensitivity band of the 
detector, M > 12OM0. For lower masses, longer waveforms 
are required, and in general can be produced by connect- 
ing post-Newtonian (PN) and numerical -relativity (NR) wave- 
forms 1 36, 65, 67, 117, 118|. The accuracy estimates given 
in this paper tell us nothing about the accuracy of such longer 
"hybrid" waveforms, because we cannot properly quantify the 
accuracy of the PN approximants. We defer the discussion 
of the accuracy of hybrid waveforms, and the implications 
for the necessary length of numerical waveforms, to separate 
work |fTT9| . 

For now we consider the fidelity of PN results to full gen- 
eral relativity only in the regime where we also have NR re- 
sults, i.e., in the last orbits before merger. We compare the PN 
and NR phase disagreement over the last 8-10 GW cycles be- 
fore Mco = 0.1 for two classes of PN approximant, TaylorTl 
and TaylorT4. For nonspinning cases we find that the perfor- 
mance of both approximants does not change drastically as the 
mass ratio is increased to ^ = 4, and this means that the Tay- 
lorT4 approximant continues to provide the best agreement, 
with an accumulated phase disagreement in the q — 4 case of 
0.2 rad, or 0. 1 rad if we consider the root-mean-square average 
of the phase disagreement; see Fig. [8] For spinning binaries, 
the two approximants include spin terms up to only 2.5PN or- 



der The TaylorTl approximant nonetheless is fairly robust, 
while TaylorT4-truncated performs poorly for large spins, in 
particular large spins anti-aligned with the binary's orbital an- 
gular momentum. The full TaylorT4 approximant performs 
well for all spins larger than Xi ~ —0.75. 

Finally, we study the accuracy of the PN wave amplitude, 
and find that when the highest order amplitude corrections are 
included (3PN for nonspinning binaries, and 2PN for spinning 
cases), the amplitude error is no more than 4%. This is in con- 
trast to the quadrupole amplitude, which can over-estimate the 
true physical amplitude by up to 13% at black hole dimen- 
sionless spins of Xi = +0.85 corresponding to an increase of 
44% in detection rates. Note that precisely the cases with the 
largest SNR (spins aligned with the angular momentum) are 
also those with the largest PN amplitude errors. 
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Appendix A: Apparent-horizon and puncture estimates of the 
black-hole masses 



There are two methods that are commonly used to estimate 
the masses of black holes in puncture data, in addition to an- 
alyzing apparent horizons. The first, which is generally ap- 
plicable to all black-hole data, is to make use of the area of 
the apparent horizon, A. The black hole's "irreducible mass" 
Mirr is given by Mirr — ^/aJI&k, and the total mass can be 
estimated by M120I 



M^ = Ml 



S^ 



4M,?. 



(Al) 



A second method is to make use of the asymptotic prop- 
erties of the wormhole puncture data. Each puncture rep- 
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resents an extra asymptotically flat end of the slice, and the 
ADM mass calculated at each "extra" end can be considered 
as a measure of the mass of that black hole. In the puncture- 
data construction, the momentum constraint is solved analyt- 
ically by the Bowen-York conformal extrinsic curvature, and 
the Hamiltonian constraint is solved numerically to give the 
function u in the ansatz |fT9ll . 



V/=l 



2n 



m2 
2^ 



(A2) 



where m, parametrizes the mass of the /th black hole, and r, 
is the coordinate distance to the rth black hole. The resulting 
data represent two black holes on a three-sheeted topology. 
One sheet contains two black holes, and represents the physi- 
cal space that we want to describe. Each black hole has an ex- 
tra sheet associated with it, which extends to an extra asymp- 
totically flat end, and in the puncture construction those ends 
are compactified to points, or "punctures". 

To calculate this mass, we require only the value of the 
function u at the puncture. The mass is then given by 



Mi = nii (1+Ui 



2D J 



(A3) 



where D is the coordinate distance between the two punctures. 
A derivation of this expression is given in 1 19 1. The two mea- 
sures of the mass that we have discussed are shown to agree 
within numerical uncertainty in the case of nonspinning black 
holes in I.121J . Since the ADM mass at the puncture can be 
easily calculated directly from the initial data with high pre- 
cision, it has become a standard tool in assessing the mass of 
black holes in puncture data. 



However, as discussed in (ST], this is only a reasonable 
measure of the black-hole mass for nonspinning black holes. 
A heuristic explanation for this effect is that the fall-off of the 
extrinsic curvature for a boosted Bowen-York black hole is 
far faster towards the extra asymptotically flat ends as it is to- 
wards the "physical" end, and so the extra sheets of the topol- 
ogy contain far less junk radiation than the physical sheet, and 
the ADM mass of each of those sheets is not contaminated 
by very much junk radiation. In the spinning case, however, 
the fall-off on the second sheet is the same as on the physical 
sheet, and so the extra sheets each contain roughly the same 
junk radiation as the physical space, and only for low spins 
will the ADM mass at the puncture be a good measure of the 
black-hole mass. 



As an illustration of this effect, the values of the black-hole 
mass as given by the two methods are shown in Table |V1| 



For the simulations presented in this paper, the results were 
first produced using the puncture-mass estimates. They were 



then rescaled according to the results in Tab. VI A rescaling 
of mass will have an overall effect on the time-scale of the 
simulations, but we found that even in the highest spin case 
the effect was negligible. This is most important in the com- 



parison with PN approximants in Sec. Vll where the PN and 
NR results are compared assuming the same mass scale. But 
we find that the phase disagreement between the NR and NR 
results is much larger than the error due to using the incorrect 
black-hole mass, and does not noticeably alter the results in, 
for example. Fig. IT] 
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